The X-Philes (2nd Revision)/The X-Philes Number 1 (1995).iso/xphiles/hp48hor2/hamilton.doc

Text File | 1995-03-31 | 10KB | 213 lines

(Comp.sys.hp48)
Item: 693 by [e-mail]. [Charles Patton]
Subj: Hamiltonian Method for Geodesics
Date: Wed Feb 19 1992

@ Version 1.0 @

_Introduction_

     The programs contained in this document comprise a collection
of simple-minded utilities for use in computing approximate orbits
of Hamiltonian dynamical systems in three space variables.

     These utilites are the outcome of messing around with ways of
computing geodesics on implicitly-defined surfaces embedded in Euclidean
space. There is nothing new about the method involved, rather, the nicest
thing is how easy it is to "follow your nose" when working in an
environment like that of the '48.

     All suggestions, additions, corrections, insights, commentary,
and criticisms are welcome.

_NO WARRANTY_
    This work is provided on an "as is" basis. Hewlett-Packard Co.
provides no warranty whatsoever, either express or implied regarding
the work, including warranties with respect to its merchantability or
fitness for any purpose whatsoever.

_Copyright_
   Copyright (C), 1991, Hewlett-Packard Co. Permission to copy all or
part of this work is granted provided that the copies are not made or
distributed for resale (excepting nominal copying fees) and the
 _NO WARRANTY_ and this copyright notice are included verbatim. Other
permissions can be arranged by contacting the author.

_Correspondence_

   Correspondence on HamiltonGeodesy should be sent to Charles Patton
at one of the following addresses:

snail-mail:
          M.S. 5U-L9
          Hewlett-Packard Co.
          1000 N.E. Circle Blvd.
          Corvallis, OR 97330

e-mail:
          [e-mail]    (for Internet hosts)
          hplabs!hpcvrs!charliep (for UUCP hosts)

_Overview_

     Recall that the Hamilton-Jacobi formulation of classical mechanics
starts with the energy function ("Hamiltonian") as a function of position
and momentum of the system in question: q,p --> H(q,p). The Hamilton
equations then say that the time derivative of the position is the
derivative of H with respect to momentum, and the time derivative of
the momentum is the negative of the derivative of H with respect to
the position: q'=dH/dp ; p'=-dH/dq. The orbits of the Hamiltonian are
the integral curves of this system of differential equations, that is,
curves q(t),p(t) which satisfy dq/dt=dH/dp and dp/dt=-dH/dq.

    Perhaps the simplest way to formulate the problem of finding
geodesics on a Riemannian space is as a Hamiltonian system. In this case,
the Hamiltonian is just the function which assigns to a position and
momentum half the squared length (according to the local metric) of the
momentum:

H(q,p)=0.5 Sum g  (q) p  p
            i,j ij     i  j

where g(q) is the metric at q. In this case, the projection onto
position space of the orbits of the Hamiltonian are the geodesics.

    If we are trying to find geodesics on a surface emdedded in
Euclidean space, we could find local coordinates for the surface,
compute the metric in local coordinates, and then use the above
formulation to find geodesics in terms of these local coordinates.
This is almost always difficult to do.

    We will take a different approach here, the key difference being
that we are assuming that the surface is given to us as the zero-set
of a function F: q --> F(q) defined on the Euclidean space. We can use
F to define a potential function on the Euclidean space with a minimum
along the desired surface and use the Euclidean metric together to define
the Hamiltonian:

          2    2          2
H(q,p)=k F(q) +0.5 ||p||

for some "large" constant, k.

    If we start with any initial q,p with F(q)=0 and p "tangent" to
the the surface, the projection onto position space of the resulting
orbit will approximate a geodesic on the desired surface. You can
think of this as describing a particle attached to the surface by a
sliding, frictionless spring. If the spring is very stiff, the
particle will follow the path of a geodesic when it is given a shove
along the surface. If the spring is not so stiff, the particle will
combine near-geodesic motion with oscillations normal to the surface.

_Organization and Descriptions of the Programs_

    The utilities are arranged in a directory with the two top-level
programs first, and the supporting utilities and variables afterwords.
The space coordinates are represented by the variables, 'x' 'y' and 'z'
(note the lower case.) The corresponding momentum coordinates are 'px'
'py' and 'pz'.

SetUpHamiltonian:
    Given the function kF, on the stack, as an expression in the space
variables, 'x' 'y' and 'z', this program creates the Hamiltonian
described above and stores it in the variable 'H' for future
reference. It then computes the components, 'xdot'...'pzdot',
corresponding to the Hamiltonian equations. It then calls on the
utility Lie4thOrder to compute a (formal) fourth order polynomial
approximation to the solution of the system of Hamilton equations,
storing the formulae for the increments as '<delta>x'...'<delta>pz'.
This takes a loooong time so go out for a coffee break after firing it
up.

After this finishes, you are ready to plot approximate geodesics on your
surface.

    You can easily customize this program to set up for an arbitrary
Hamiltonian in these variables by deleting the first line of the
program listed, below, and storing your choice of Hamiltonian (as an
expression in 'x'...'pz') in the variable 'H' before running
SetUpHamiltonian.

N.B.>Insure that 'x' 'y' 'z' 'px' 'py' and 'pz' are formal (no variables
N.B.>with those names along the current path.)

N.B.> Examine the MODES menu to be sure that the
N.B.> the machine is in SYMBOLIC EVALUATION mode ( SYM button.)

RunHamiltonian:
    This program takes initial values for x y z px py and pz on the
stack and computes successive points (approximately) on the orbit of
the Hamiltonian. It graphs these by lines connecting successive (x,z)
pairs (i.e. projection of the orbit on the x-z plane) using the
current PPAR values for the horizontal and vertical scaling of the
screen. It requires that you have previously run SetUpHamiltonian in
the current directory.

    As part of its initialization, it takes the current value of the
stepsize (stored in the variable <delta>) and pre-computes the higher
order stepsizes (<delta>^2/2, <delta>^3/6, and <delta>^4/24) storing
these in <delta>2 <delta>3 and <delta>4. As explained below, the
default value of <delta> is 0.1, but you can use any value you wish as
long as you store it before running this program.

    You should exit this program by pressing the [ENTER] key. When
this is done, the program will clean up and leave the last computed
values of x...pz on the stack ready for continuation of the orbit, if
desired.

HamiltonStep:
    This simple utility takes the values off the stack, stores them
in 'x'...'pz', computes the next step and returns the values to the
stack.

Lie4thOrder:
    This routine constitutes the heart of this collection of
utilities.  The idea is that applying the Lie derivative operator

     d      d      d       d      d      d
L:=x'-  + y'-  + z'-  + px'- + py'- + pz'-
     dx     dy     dz     dpx    dpy    dpz

to any function, f, corresponds to taking the time derivative of f
along the orbits of the Hamiltonian. Thus, the fourth-order Taylor
series approximation to the time evolution of f is

    f+(Lf)t+(LLf)t^2/2+(LLLf)t^3/6+(LLLLf)t^4/24

    Using this on each of the coordinate functions, x...pz, thus gives
an approximation to the orbits themselves. This program computes this
approximation for each of the coordinate functions using <delta> as
the time-step, t. This is equivalent (in the limit of small time-step)
to the fourth-order Runge-Kutte method of approximating solutions to
ordinary differential equations, but has the advantage of being much
more comprehensible. It also allows you to compute the approximate
time evolution of any smooth function, not just the coordinate
functions, in exactly the same manner.

   The resulting approximation is subject to the same kinds of
instabilities as the R-K method. In our case, this shows up when the
approximation overshoots into a region where the added potential is
steep. This causes the next step to overshoot even more in the other
direction, quickly leading to huge values for the position and
momentum coordinates. For this reason, you should leave <delta> at its
default value, and if instabilites become apparent, choose a smaller
value for the initial velocity. The point here being that while the
space-projection of the orbits of the geodesic equation don't depend
on the initial velocity, the orbits of our approximation do. They are
only asymptotic to the geodesics in the limit of small initial
velocity. Moreover, no choice of potential will fix that. You can
think of it like a motorcycle driving in a groove. At low velocity, it
will follow the groove nicely, but at high velocity, it will "bank"
its turns on the wall of the groove, hence the possibility for
overshoot becomes more serious.

LieDer:
    This utility takes a function, given as an expression in x...pz,
and computes its Lie derivative according to the current values of
'xdot'...'pzdot' which correspond to x'...pz'.

The Other Variables:
    <delta>...<delta>4 and 'H' were explained previously. The others
are placeholders for values computed during the operations of the
main programs.